### MC simulations to examine legislative influence
### Simulations to show that Multiple Imputation may help alleviate the some of the measurement error issues


rm(list=ls())
library(simex)
library(arm)
library(msm)

nissues<-150
nsim<-1000
merror<-8
ncouncil<-7


difftrueobs<-matrix(NA,nrow=nsim,ncol=4)
colnames(difftrueobs)<-c("EPmean","Councilmean","EPsd","Councilsd")

disttrueout<-matrix(NA,nrow=nsim,ncol=2)
colnames(disttrueout)<-c("EP","Council")

mc.reg<-matrix(NA,nrow=nsim,ncol=6)
colnames(mc.reg)<-c("eptrue","ctrue","epobs","c.obs.mean","simex.ep","simex.c")
	
for (j in 1:nsim) {
		## Generate True Positions and Outcomes
		## Council True Position is mean/median of Council
		council<-matrix(round(runif(ncouncil*nissues,min=0,max=100),0),nrow=nissues, byrow=T)
		#True council position (mean)
		mc<-apply(council,1,mean)	
		# EP True Position is drawn from uniform dist somewhat more integrationist than mem states
		ep<-round(runif(nissues,40,100),0)

		# Outcome is equidist compromise plus stochastic error
		out<-0.5*ep+0.5*mc+rnorm(nissues,0,5)

		## Generate observed positions and outcome
		## Council observed position is mean of observed Council positions
		# observed positions of 7 council members are drawn from a random normal dist with the mean equal to true position and sd= 10
		c.obs<-council+matrix(rnorm(ncouncil*nissues,0,merror),nrow=nissues,byrow=T)
			
		## EP observed drawn from normal dist centered on EP true position, with sd =7 (same as council)
		ep.obs<-rnorm(nissues,ep,merror)
		
		difftrueobs[j,1]<-mean(abs(ep-ep.obs))
		difftrueobs[j,2]<-mean(abs(mc-(apply(c.obs,1,mean))))

		difftrueobs[j,3]<-sd(abs(ep-ep.obs))
		difftrueobs[j,4]<-sd(abs(mc-(apply(c.obs,1,mean))))

		disttrueout[j,1]<-mean(abs((0.5*ep+0.5*mc)-ep.obs))
		disttrueout[j,2]<-mean(abs((0.5*ep+0.5*mc)-(apply(c.obs,1,mean))))


		## Observed outcome drawn from normal dist centered on true outcome
		out.obs<-rnorm(nissues,out,merror)

		#rescale all council mems, EP, and outcome to place observed obs on 0-100 scale
		absmax <- max(apply(cbind(c.obs,ep.obs,out.obs),1,max))
		absmin<- min(apply(cbind(c.obs,ep.obs,out.obs),1,min))
		range<- absmax-absmin
		c.obs <- round(((c.obs+ abs(absmin))/range)*100,0)
		ep.obs <- round(((ep.obs+ abs(absmin))/range)*100,0)
		out.obs <- round(((out.obs+ abs(absmin))/range)*100,0)

		avc.obs<-	apply(c.obs,1,mean)	# council as mean of observed
		mc.obs<-	apply(c.obs,1,median) # council as median of observed
	
		true<-lm(out ~ ep+mc-1)
		mc.reg[j,1]<-true$coefficients[1]
		mc.reg[j,2]<-true$coefficients[2]

		obs<-lm(out.obs~ep.obs+avc.obs-1)
		mc.reg[j,3]<-obs$coefficients[1]
		mc.reg[j,4]<-obs$coefficients[2]
		
		simex1<-simex(model=obs,SIMEXvariable="ep.obs",measurement.error=6.39,asymptotic=F)
		mc.reg[j,5]<-simex1$coefficients[1]
		mc.reg[j,6]<-simex1$coefficients[2]

}

save(mc.reg,file="~/Dropbox/ThomsonDEU/simout_25Feb13.RData")

colMeans(mc.reg)
			
library(sfsmisc)

boxplot.matrix(mc.reg, names=c("EP","Council","EP","Council","EP","Council"),ylab="Coefficients",ylim=c(0.25,.9))
text(x = 1.5, y = 0.8, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=1.5, y=.85, "Model 1 \nTrue")
text(x = 3.5, y =  0.8, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=3.5, y=.85, "Model 2 \nObserved (Mean)")
text(x = 5.5, y = 0.8, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=5.5, y=0.85, "Model 4 \nSimex")


outmean1<-mean(ifelse(mc.reg[,2]> mc.reg[,1],1,0)) 
outmean2<-mean(ifelse(mc.reg[,4]> mc.reg[,3],1,0)) 
outmean3<-mean(ifelse(mc.reg[,6]> mc.reg[,5],1,0))


colMeans(difftrueobs)
colMeans(disttrueout)





